
package org.zhongyin.Regression;
import java.io.IOException;
import java.text.ParseException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

import jxl.read.biff.BiffException;

import org.joda.time.DateTime;
import org.zhongyin.Util.Utility;
import org.zhongyin.internal.core.ReadExcel;

import Jama.Matrix;


public class Regression {
    
    private static final double TUNING_CONSTANT = 4.685d;
    
    private Matrix coefficient;
    
    private double leastSquares;
    
    private double standardDeviation;
    
    private Matrix weight;
    
    //D = sqrt(eps(class(double))) from STATROBUSTFIT_Revised.m under matlab
    private static final double D = 1.490116119384766e-08;
    
    private static final int iterLimit = 50;
    
    public Regression(){
            coefficient = null;
            
            leastSquares = 0d;
            
            standardDeviation = 0d;
            
            weight = null;
            
    }
    
    
    
    /**
     * For Linear Regression
     * b = (XT*X)-1 * XT * Y
     * @param array_X       X
     * @param array_Y       Y
     * @return      coefficient (b)
     */
    public Matrix getCoefficient(double[][] array_X, double[] array_Y){
            
            if(array_X == null || array_Y == null){
                    throw new NullPointerException();
            }
            Matrix X = new Matrix(array_X);
            Matrix Y = new Matrix(array_Y, array_Y.length);
            
            Matrix X_Transpose = X.transpose();
            // XT*X
            Matrix X_Transpose_Times_X = X_Transpose.times(X);
            //(XT*X)-1 Matrix inverse
            Matrix square_inverse = X_Transpose_Times_X.inverse();
            //(XT*X)-1 * XT
            Matrix inverseTimesTranspose = square_inverse.times(X_Transpose);
            
            Matrix result = inverseTimesTranspose.times(Y);
            
            
            coefficient = result;
            
            return result;
            
    }
    
    
    // for linear regression
    /**
     * @param array_Y       Y
     * @param array_X       X
     * @param array_b       coefficient
     * @return      leastsquares (e2)
     */
    public double getLeastSquares(double[] array_Y, double[][] array_X, double[][] array_b){
            if(array_Y == null || array_X == null || array_b == null){
                    throw new NullPointerException();
            }
            
            Matrix X = new Matrix(array_X);
            Matrix b = new Matrix(array_b);
            Matrix Xb = X.times(b);
            
            
            double[][] array_Xb = Xb.getArray();
            
            double[][] e = new double[array_Y.length][1];
            
            
            double result = 0d;
            
            
            for(int i = 0; i < e.length; i++){
                    
                    e[i][0] = array_Y[i] - array_Xb[i][0];
                    result += e[i][0] * e[i][0];
            }
            
            
            
            leastSquares = result;
            
            return result;
            
            
    }
    
    
    /**
     * ri = (yi - ybari) / std
     * @param array_Y       Y
     * @param array_X       X
     * @param array_b       coefficient
     * @param sDeviation 
     * @return      Standardized residuals
     */
    public Matrix getStdResidual(double[] array_Y, double[][] array_X, double[][] array_b, double sDeviation){
            if(array_Y == null || array_X == null || array_b == null){
                    throw new NullPointerException();
            }
            
            Matrix X = new Matrix(array_X);
            Matrix b = new Matrix(array_b);
            Matrix Xb = X.times(b);
            
            double[][] y_bar = Xb.getArray();
            
            double[] r = new double[array_Y.length];
            
            for(int i = 0; i < r.length; i++){
                    r[i] = (array_Y[i] - y_bar[i][0])/sDeviation;
            }
            
            Matrix stdResidual = new Matrix(r, array_Y.length);
            
            return stdResidual;
    }
    
    
    
    /**
     * 
     * Compute the adjusted residuals
     * @param array_Y Y
     * @param array_X X
     * @param array_b coefficient
     * @param array_h leverage
     * @return
     */
    public double[] getAdjustedResiduals(double[] array_Y, double[][] array_X, double[][] array_b, Matrix array_h){
            if(array_Y == null || array_X == null || array_b == null){
                    throw new NullPointerException();
            }
            
            Matrix X = new Matrix(array_X);
            Matrix b = new Matrix(array_b);
            Matrix Xb = X.times(b);
            
            double[][] h = array_h.getArray();
            
            double[][] y_bar = Xb.getArray();
            
            double[] r = new double[array_Y.length];
            
            for(int i = 0; i < r.length; i++){
//                  System.out.print(array_Y[i] - y_bar[i][0]+ ", ");
                    r[i] = (array_Y[i] - y_bar[i][0])/Math.sqrt(1-h[i][i]);
//                  r[i] = 1/Math.sqrt(1-h[i][i]);
            }
            
//          System.out.println();
            
            return r;
    
    }
    
    /**
     * s is the robust variance given by MAD/0.6745 where MAD 
     * is the median absolute deviation of the residuals
     * 
     * @param array_Radj  Radj
     * @param array_X X
     * @return
     */
    public double getS(double[] array_Radj, double[][] array_X){
            if(array_Radj == null){
                    return 0;
            }
            
            
             double[] array_Radj_Clone = array_Radj.clone();
            
             double median = 0d;
             
         int length = array_X[0].length; 
         
//       System.out.println ("length : " + length);
         for(int i = 0; i < array_Radj.length; i++){
             array_Radj_Clone[i] = Math.abs(array_Radj[i]);
         }
            
             Arrays.sort(array_Radj_Clone);
             
             
             
             array_Radj_Clone = Arrays.copyOfRange(array_Radj_Clone, length-1, array_Radj.length);
             
//           System.out.println("Sorted Radj: " + Arrays.toString(array_Radj_Clone));
             
             if(array_Radj.length % 2 == 0){
                     
                     
                     median = (array_Radj_Clone[array_Radj_Clone.length/2] + array_Radj_Clone[array_Radj_Clone.length/2 - 1])/2;
             }
             else{
                     median = array_Radj_Clone[array_Radj_Clone.length/2];
             }
             
//           median = array_Radj[array_Radj.length/2];
             
//           System.out.println("Media :" + median);
             double MAD = median;
             
             double S = MAD / 0.6745;
             
             return S;
            
            
    }
    
    /**
     * 
     * @param stdResidual ri
     * @return      W[][] with Wii(ri)
     */
    public Matrix getOriginalWeight(Matrix stdResidual){
            if(stdResidual == null){
                    throw new NullPointerException();
            }
            
            double[][] r = stdResidual.getArray();
            
            double[][] W = new double[r.length][r.length];
            
            for(int i = 0; i < r.length; i++){
                    
                    if(Math.abs(r[i][0]) <= TUNING_CONSTANT){
                            W[i][i] = (1 - (r[i][0] * r[i][0]) / (TUNING_CONSTANT * TUNING_CONSTANT))*(1 - (r[i][0] * r[i][0]) / (TUNING_CONSTANT * TUNING_CONSTANT));
                    }
                    else{
                            W[i][i] = 0d;
                    }
            }
            
            
            
            Matrix Weight = new Matrix(W);
            
            weight = Weight;
            
            return Weight;
            
            
    }
    
    public Matrix getRobustWeight(Matrix U){
            if(U == null){
                    throw new NullPointerException();
            }
            
            double[][] u = U.getArray();
            
            double[][] W = new double[u.length][u.length];
            
            for(int i = 0; i < W.length; i++){
                    
                    if(Math.abs(u[i][0]) <= 1){
                            W[i][i] = (1 - (u[i][0] * u[i][0]))*(1 - (u[i][0] * u[i][0]));
                    }
                    else{
                            W[i][i] = 0d;
                    }
            }
            
            
            
            Matrix Weight = new Matrix(W);
            
            weight = Weight;
            
            return Weight;
            
            
    }
    
    
    /**
     * 
     * The quantities hi are the leverages
     * , the diagonal elements of X(Xt*X)-1 * Xt
     * @param array_X X
     * @return array_h
     */
    public Matrix getLeverage(double[][] array_X){
            if(array_X == null){
                    throw new NullPointerException();
            }
            
            Matrix X = new Matrix(array_X);
            
            // h = X(XTX)-1XT
            Matrix result = X.times((X.transpose().times(X)).inverse()).times(X.transpose());
            
            return result;
    }
    
    
    /**
     * @param leastSquares e2
     * @return e
     */
    public double getStandardDeviation(double leastSquares){
            
            double result = Math.sqrt(leastSquares);
            
            standardDeviation = result;
            
            return result;
            
    }
    
    
    
    /**
     * For Weighted Linear Regression
     * b = (XT*WX)-1 * XT * W *Y
     * @param array_X       X
     * @param array_Y       Y
     * @return      coefficient (b)
     */
    public Matrix getWeightedCoefficient(double[] array_Y, double[][] array_X, Matrix weight){
            if(array_X == null || array_Y == null || weight == null){
                    throw new NullPointerException();
            }
            Matrix X = new Matrix(array_X);
            Matrix Y = new Matrix(array_Y, array_Y.length);
            
            Matrix X_Transpose = X.transpose();
            
            //XT*W
            Matrix X_Transpose_Times_W = X_Transpose.times(weight);
            // XT*W*X
            Matrix X_Transpose_Times_W_Times_X = X_Transpose_Times_W.times(X);
            //(XT*WX)-1 Matrix inverse
            Matrix square_inverse = X_Transpose_Times_W_Times_X.inverse();
            //(XT*WX)-1 * XT
            Matrix inverseTimesTranspose = square_inverse.times(X_Transpose);
            //(XT*WX)-1 * XTW
            Matrix timesWeight = inverseTimesTranspose.times(weight);
            //(XT*WX)-1 * XTWY
            Matrix result = timesWeight.times(Y);
            
            
            coefficient = result;
            
            return result;
    }
    
    
    
    
    /**
     * For Weighted Linear Regression
     * @param array_Y       Y
     * @param array_X       X
     * @param array_b       coefficient
     * @return      leastsquares (e2)
     */
    public double getWeightedLeastSquares(double[] array_Y, double[][] array_X, double[][] array_b, Matrix weight){
            if(array_Y == null || array_X == null || array_b == null || weight == null){
                    throw new NullPointerException();
            }
            
            Matrix X = new Matrix(array_X);
            Matrix b = new Matrix(array_b);
            Matrix Xb = X.times(b);
            
            
            double[][] array_Xb = Xb.getArray();
            
            double[][] e = new double[array_Y.length][1];
            
            
            double result = 0d;
            
            
            for(int i = 0; i < e.length; i++){
                    
                    e[i][0] = array_Y[i] - array_Xb[i][0];
                    //Wii*(Yi-Xib);
                    result += weight.getArray()[i][i] * (e[i][0] * e[i][0]);
            }
            
            
            
            leastSquares = result;
            
            return result;
            
            
    }
    
    
    
    public Matrix getU(double[] array_Radj, double S){
            if(array_Radj == null){
                    return null;
            }
            
            double[] u = new double[array_Radj.length];
            
            
            
            
            for(int i = 0; i < array_Radj.length; i++){
                    u[i] = array_Radj[i]/(TUNING_CONSTANT*S);
            }
            
            Matrix U = new Matrix(u, u.length);
            
            return U;
            
            
    }
    /**
     * Iteratively get the weight by Irls
     * @param array_Y Y
     * @param array_X X
     * @param weight array of weight
     * @return get more accurate weight
     */
    public Matrix getIterWeight(double[] array_Y, double[][] array_X, Matrix weight){
            
//          Regression regression = new Regression();
            
            int length = weight.getArray().length;
            
            int iter = 0;
            
//          //Generate an identity matrix
//          weight = Matrix.identity(length, length);
            
            double[] weight_array = new double[length];
            
            for(int i = 0 ; i < length; i++){
                    weight_array[i] = weight.get(i, i);
            }
            
//          System.out.println(iter + " weight: " + Arrays.toString(weight_array));
            // First coefficient (b)
            Matrix array_b = getWeightedCoefficient(array_Y, array_X, weight);
//          System.out.println(iter + " b: " + Arrays.deepToString(array_b.getArray()));
            // First leverage
            Matrix leverage = getLeverage(array_X);
//          System.out.println(iter + " h: " + Arrays.deepToString(leverage.getArray()));
            // First adjusted residuals
            double[] array_Radj = getAdjustedResiduals(array_Y, array_X, array_b.getArray(), leverage);
//          System.out.println(iter + " Radj: " + Arrays.toString(array_Radj));
            // First S
            double S = getS(array_Radj, array_X); 
//          System.out.println(iter + " S: " + S);
            // First U
            Matrix U = getU(array_Radj, S);
//          System.out.println(iter + " U: " + Arrays.deepToString(U.getArray()));
            
            
            
            
            
//          System.out.println();
            
            
            
            while(any(array_b.getArray()) && iter < iterLimit){
                    
                    iter++;
                    
                    weight = getRobustWeight(U);
                    
                    for(int i = 0 ; i < length; i++){
                            weight_array[i] = weight.get(i, i);
                    }
                    
//                  System.out.println(iter + " weight: " + Arrays.toString(weight_array));
                    
                    array_b = getWeightedCoefficient(array_Y, array_X, weight);
//                  System.out.println(iter + " b: " + Arrays.deepToString(array_b.getArray()));
                    // First leverage
                    leverage = getLeverage(array_X);
//                  System.out.println(iter + " h: " + Arrays.deepToString(leverage.getArray()));
                    // First adjusted residuals
                    array_Radj = getAdjustedResiduals(array_Y, array_X, array_b.getArray(), leverage);
//                  System.out.println(iter + " Radj: " + Arrays.toString(array_Radj));
                    // First S
                    S = getS(array_Radj, array_X); 
//                  System.out.println(iter + " S: " + S);
                    // First U
                    U = getU(array_Radj, S);
//                  System.out.println(iter + " U: " + Arrays.deepToString(U.getArray()));
                    
                    
                    
                    
                    
//                  System.out.println();
                    
                    
            }
            
            
            double[][] diagonal_weight = new double[weight_array.length][weight_array.length];
            
            for(int i = 0 ; i < weight_array.length; i++){
                    diagonal_weight[i][i] = weight_array[i];
            }
            Matrix result = new Matrix(diagonal_weight);
            
            return result;
    }
    
    /**
     * 
     * True if any element of a vector is
     * logical 1 (TRUE).
     * @param array
     * @return
     */
    public boolean any(double[][] array){
            if(array == null){
                    return false;
            }
            
            for(int i = 0; i < array.length; i++){
                    if(Math.abs(array[i][0])  <= D * Math.max(Math.abs(array[i][0]), 0)){
                            return false;
                    }
            }
            
            return true;
    }
    
    
    /**
     * get forgetting factor based on lumda from 0.970 to  1.000
     * @param weight
     * @return a list of forgetting factor
     */
    public List<double[][]> getWeightPrimePrime(Matrix weight){
            if(weight == null){
                    throw new NullPointerException();
            }
            
            double k = 0.971d;
            
            List<double[][]> list = new ArrayList<double[][]>();
            
            int length = weight.getArray().length;
            
            
            while( k <= 0.971){
                    Matrix lumda = Matrix.identity(length, length);
                    double[][] array_lumda = lumda.getArray();
                    
                    System.out.print("w: ");
                    for(int i = 0; i < length; i++){
                            array_lumda[i][i] = lumda.get(i, i) * Math.pow(k, length-1 - i) * weight.get(i, 0);
                            System.out.print(array_lumda[i][i]+ ", ");
                    }
                    

                    list.add(array_lumda);
                    k += 0.001;
            }
            
            return list;
    }
    
    public double getSecondStageWeightedLeastSquare(double[] array_Y, double[][] array_X, List<double[][]> weightPrimePrime){
            if(array_Y == null || array_X == null || weightPrimePrime == null){
                    throw new NullPointerException();
            }
            
            double minE = Double.MAX_VALUE;
            double minLeastSquare = Double.MAX_VALUE;
            
            int counter = 0;
            for(double[][] array_weight: weightPrimePrime){
                    
                    if(counter == 29){
                            System.out.println(counter + "weight: " + Arrays.deepToString(array_weight));
                    }
                    
//                  System.out.println(counter + "weight: " + Arrays.deepToString(array_weight));
                    
                    Matrix weight = new Matrix(array_weight);
                    Matrix coefficient = getWeightedCoefficient(array_Y, array_X, weight);
                    
                    double[][] array_b = coefficient.getArray();
                    
                    System.out.println(counter+ "b:" + Arrays.deepToString(array_b));
                    
                    double weightedLeastSquare = getWeightedLeastSquares(array_Y, array_X, array_b, weight);
                    
                    double e = Math.sqrt(weightedLeastSquare/array_Y.length);
                    
                    System.out.println(counter+ "e2: " + weightedLeastSquare);
                    
                    System.out.println(counter+ "e: " + e);
                    
                    minE = Math.min(minE, e);
                    
                    minLeastSquare = Math.min(minLeastSquare, weightedLeastSquare);
                    
                    counter++;
                    
            }
            
            return minE;
    }
    
    
    public List<Matrix> getFFCofficient(double[] array_Y, double[][] array_X, List<double[][]> weightPrimePrime){
            if(array_Y == null || array_X == null || weightPrimePrime == null){
                    throw new NullPointerException();
            }
            
            double lumbda = 0.971d;
            
            
            
            List<Matrix> result = new ArrayList<Matrix>(); 
            for(double[][] array_weight: weightPrimePrime){
                    
                    
                    

                    Matrix weight = new Matrix(array_weight);
                    Matrix coefficient = getWeightedCoefficient(array_Y, array_X, weight);
                    
                    System.out.println();
                    
                    System.out.println(lumbda + " b: " + Arrays.deepToString(coefficient.getArray()));
                    
                    result.add(coefficient);
                    
                    
                    lumbda += 0.001;
                    
            }
            
            return result;
            
            
    }
	
	public static void main(String args[]) throws IOException, ParseException, BiffException{
		double[][] array_X = {{1,2},{1,4},{1,7}};
		double[] array_Y = {1,3,5};
		double[] array_Weight = {1,1,1};
		
		int length = 3;
//		
//		ReadExcel Original = new ReadExcel();
//		Original.setInputFile("c:/CoolingLoadProject/All result1.xls");
//		
//	    double[][] original_data = Original.readT24_168(1);
//	    
//	    int length_origianl = original_data.length;
//	    double[] array_Y = new double[214];
//	    double[][] array_X = new double[214][17];
//	    
//	    
//	    
//	    for(int i = 5, k = 0; i < 5140&&k<214; i=i+24, k++){
//	    	array_Y[k] = original_data[i][0];
//	    	for(int j = 0; j < 11; j++){	
//	    		array_X[k][j] = original_data[i][j+1];
//	    	}
//	    }
//	    
//	    for(int i = 5, k = 0; i < 5140&&k<214; i=i+24, k++){
//	    		//t*t
//	    		array_X[k][11] = original_data[i][6] * original_data[i][6];
//	    		//t24*t24
//	    		array_X[k][12] = original_data[i][7] * original_data[i][7];
//	    		//t168*t168
//	    		array_X[k][13] = original_data[i][8] * original_data[i][8];
//	    		//t * RHt
//	    		array_X[k][14] = original_data[i][6] * original_data[i][9];
//	    		//t24 *RHt24
//	    		array_X[k][15] = original_data[i][7] * original_data[i][10];
//	    		//t168 * RHt168
//	    		array_X[k][16] = original_data[i][8] * original_data[i][11];
//	    	
//	    }
	    
	    
	    
//	    Matrix Y = new Matrix(array_Y, array_Y.length);
	    
//	    double[] test_Y = Y.getRowPackedCopy();
	    
//	    System.out.println(Arrays.toString(test_Y));
	    
	    
	    
	    Utility.printArray(array_X, 0, array_X.length);
	    
	    System.out.println(Arrays.toString(array_Y));

	   
//	    
//	    int length = array_Y.length;
	    
	    Matrix weight = Matrix.identity(length, length);
		
		Regression regression = new Regression();
		
		weight = regression.getIterWeight(array_Y, array_X, weight);
		
		Utility.printArraySingleLine(weight.getArray());
		
		Matrix cofficient = regression.getWeightedCoefficient(array_Y, array_X, weight);
		
		System.out.println("ARXModel(24, 168) b: " + Arrays.deepToString(cofficient.getArray()));
				
		//		List<double[][]> weightPrimePrime = new ArrayList<double[][]>();
				
		//		weightPrimePrime = regression.getWeightPrimePrime(weight);
				
		//		double minLeastSquare = regression.getSecondStageWeightedLeastSquare(array_Y, array_X, weightPrimePrime);
				
		//		System.out.println("Min Least Square: " + minLeastSquare);
				
		//		regression.getFFCofficient(array_Y, array_X, weightPrimePrime);
		
		/**
		ARXFirstWDProvider test = new ARXFirstWDProvider();
		test.setHistdata();
		for(int i = 1; i <= 6; i++){
			
				test.setARXHistmodel(i, 0);
		//		
//				Utility.printArraySingleLine(test.getARXmodel_X());
//			
//				System.out.println(Arrays.toString(test.getARXmodel_Y()));
				
				double[][] array_X = test.getARXmodel_X();
				
				double[] array_Y = test.getARXmodel_Y();
				
				int length = array_Y.length;
			    
			    Matrix weight = Matrix.identity(length, length);
				
				Regression regression = new Regression();
				
				weight = regression.getIterWeight(array_Y, array_X, weight);
				
//				Utility.printArraySingleLine(weight.getArray());
				
				Matrix cofficient = regression.getWeightedCoefficient(array_Y, array_X, weight);
				
				System.out.println("ARXModel("+i+"," + (i+1) + ", 24) b: " + Arrays.deepToString(cofficient.getArray()));
				}
		
		ARXSecondWDProvider second = new ARXSecondWDProvider();
		
		
		second.setHistdata();
		second.setARXHistmodel(0);
		//		
//		Utility.printArraySingleLine(second.getARXmodel_X());
//	
//		System.out.println(Arrays.toString(second.getARXmodel_Y()));
		
		double[][] array_X = second.getARXmodel_X();
		
		double[] array_Y = second.getARXmodel_Y();
		
		int length = array_Y.length;
	    
	    Matrix weight = Matrix.identity(length, length);
		
		Regression regression = new Regression();
		
		weight = regression.getIterWeight(array_Y, array_X, weight);
		
//				Utility.printArraySingleLine(weight.getArray());
		
		Matrix cofficient = regression.getWeightedCoefficient(array_Y, array_X, weight);
		
		System.out.println("ARXModel(24, 168) b: " + Arrays.deepToString(cofficient.getArray()));
				
		//		List<double[][]> weightPrimePrime = new ArrayList<double[][]>();
				
		//		weightPrimePrime = regression.getWeightPrimePrime(weight);
				
		//		double minLeastSquare = regression.getSecondStageWeightedLeastSquare(array_Y, array_X, weightPrimePrime);
				
		//		System.out.println("Min Least Square: " + minLeastSquare);
				
		//		regression.getFFCofficient(array_Y, array_X, weightPrimePrime);
		 */
		
	}
	
	
	
}
